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SUMMARY 


A frequency domain maximum likelihood method is developed for the estima- 
tion of airplane stability and control parameters from measured data. The 
model of an airplane is represented by a discrete-type steady-state Kalman filter 
with time variables replaced by their Fourier series expansions. The likeli- 
hood function of innovations is formulated, and by its maximization with 
respect to unknown parameters the estimation algorithm is obtained. This algo- 
rithm is then simplified to the output error estimation method with the data in 
the form of transformed time histories, frequency response curves, or spectral 
and cross-spectral densities. The development is followed by a discussion on 
the equivalence of the cost function in the time and frequency domains, and on 
advantages and disadvantages of the frequency domain approach. The algorithm 
developed is applied in four examples to the estimation of longitudinal param- 
eters of a general aviation airplane using computer-generated and measured data 
in turbulent and still air. The cost functions in the time and frequency 
domains are shown to be equivalent; therefore, both approaches are complemen- 
tary and not contradictory. Despite some computational advantages of parameter 
estimation in the frequency domain, this approach is limited to linear equa- 
tions of motion with constant coefficients. 


INTRODUCTION 

The early approaches to the extraction of airplane stability and control 
parameters from flight data were based on simple semigraphical or analytical 
methods. Some of these methods used measured frequency response curves which 
provided good insight into the physics of the system and reduced data process- 
ing to the use of simple algebra. One of the first attempts to analyze mea- 
sured data in the frequency domain for obtaining the characteristics of the 
short-period longitudinal motion of an airplane was made in reference 1 . In 
reference 2 the same characteristics were estimated either by fitting the mea- 
sured frequency response curves or by substituting the measured data in the 
transfer function equation and minimizing the resulting error. In both cases 
the least-squares technique was applied. The same technique was used for the 
direct estimation of the longitudinal and lateral aerodynamic parameters in 
references 3 and 4, respectively. 

The regression with complex variables was developed in reference 5 and 
applied to the estimation of airplane transfer function coefficients from mea- 
sured frequency response curves. A more general formulation of the regression 
in the frequency domain was introduced in reference 6 and extended to the maxi- 
mum likelihood method in reference 7. In both cases the procedure was used for 
the design of an optimal input for system identification rather than for param- 
eter estimation. 

With the availability of modern digital computers, the frequency domain 
for airplane parameter estimation was almost forgotten and the measured data 



have been mostly analyzed in the time domain. However, some further research 
and applications in this area have appeared. New frequency domain methods for 
system identification based on the equation-error formulation were introduced 
in reference 8. Frequency domain data were used for the extraction of param- 
eters of an elastic airplane in reference 9, of parameters of an airplane with 
nonsteady aerodynamics in references 10 and 11, and of flying qualities crite- 
ria in reference 12. 

The material contained in this report is an extension of the research ini- 
tiated in reference 5 and continued in references 6 and 7. Also included in 
this report are some of the developments and results from references 13 and 14, 
respectively. The purpose of this report is to present a rigorous development 
of an algorithm for the maximum likelihood estimation of airplane parameters in 
the frequency domain. The report also briefly points out the relationships 
between the estimation in the time and frequency domains, and the advantages 
and disadvantages of the frequency domain approach, mainly in terms of appli- 
cability, computing complexity, and accuracy of final results. The development 
starts with the formulation of a steady-state Kalman filter for a linear dynam- 
ical system. Before the log-likelihood function of the innovations is formu- 
lated, the basic properties of a complex random number and random sequence are 
presented. The log-likelihood function is minimized by using the modified 
Newton-Raphson technique. The maximum likelihood algorithm is then simplified 
by neglecting external disturbances to the airplane. Following the discussion, 
four examples are presented. They deal with the simplified longitudinal motion 
of a general aviation airplane and use both computer-generated and real-flight 
data. 


SYMBOLS 

A sensitivity matrix 

a z reading of vertical accelerometer, g units 

B covariance matrix of residuals 

Cjk pitching-moment coefficient, M Y /qSc 

C z vertical-force coefficient, F z /qS 

c wing mean aerodynamic chord, m 

c-| ,C 2 constants in differential equation for Gauss-Markoff process 
D matrix of transformed- system equations 

e{ ) expected value 

F matrix of continuous system 

F z force along vertical body axis, N 


2 


G 


control matrix of continuous system 
G w process-noise distribution matrix of continuous system 

g acceleration due to gravity, m/sec 2 

H transformation matrix 

I identity matrix 

Iy moment of inertia about lateral body axis, kg-m 2 

J (©) log-likelihood function 

3 = (FT 


K 


K n ,12,. 
k 1,2,. . 
M 
My 


m,n 


m 

N 

P 

P{ } 

P r ] 

Q 

q 

q 

R 

R zz 


Kalman-f ilter-gain matrix 
w/ 33 elements of K matrix 

constants in equations of motion 
Fisher information matrix 
pitching moment, N-m 

transformed quantity at mth and nth interval, respectively 

mass, kg (in appendix D) 

number of data points 

covariance matrix of state variables 

probability 

probability density 

process-noise covariance matrix 

rate of pitch, rad/sec 


kinetic pressure, 


-PV 2 , 

2 


N/m 2 


measurement-noise covariance matrix 
correlation function of z 


r 


number of output variables 



y and u 


S 

Sy U 

®ZZ 


S,t,T 


T 

u 

V 


V 

w 

w g 

X 

y 

z 


z 


n 


a 


«v 

r 



e 


0 


9 

v 

P 

a 2 


wing area, m 2 
cross-spectral density of 
spectral density of z 
quantity at sth, tth, and Tth interval, respectively 
transfer-function matrix 
control vector 
true airspeed, m/sec 
measurement-noise vector 
process-noise vector 

vertical component of turbulence velocity, nv'sec 

state vector 

measurement vector 

random variable real or complex 

= exp(jnw 0 ) 

angle of attack, rad 

angle of attack measured by wind vane, rad 
control matrix of discrete system 

process-noise distribution matrix of discrete system 
elevator deflection, rad 

Kronecker delta 

arbitrary small number 

vector of unknown parameters 

pitch angle, rad 

innovation vector 

air density, kg/m 3 

variance (a is standard deviation) 
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transition matrix 


$ 


cp y 

<Pyu 

00 


phase angle of complex variable y, deg 

phase-angle characteristics relating y and u variables, deg 
angular frequency, rad/sec 


oj q = 2rr/N 

Aerodynamic derivatives (referenced to a system of body axes with the origin at 
the airplane center of gravity) ; 


% 


Cm 


^Cjn 

2V 0 

^Cm 


5e 36, 


9 Cn 


8 Cm 
= 3c 


3C, 


% = 


qc 

2V 0 




-m 


ac 

2V 0 


'Za 


3C Z 

3a~ 


J 6e 36, 


c ro a ,c m q ' 

CSq'Cmge 


defined in appendix D (eqs. (D5) to (D8) ) 


Subscripts: 

c continuous system 

E measured quantity 

g gust 

k kth element of vector or kth column of matrix 

Z Jlth element of vector or Zth row of matrix 

m vector consisting of all elements up to and including m 

0 initial value 
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Superscripts: 

T transpose matrix 

-1 inverse matrix 

estimated value 
transformed variable 

. derivative with respect to time 

* transpose complex conjugate matrix 

R real part 

I imaginary part 

Mathematical notation: 

Tr trace of matrix 

Re real part of complex number 

I I determinant 


y 

u 


A 


amplitude-ratio characteristic relating y and u variables 
increment 


ESTIMATION ALGORITHM 

For the development of the estimation algorithm it is necessary first to 
postulate the model of an airplane and then to transform this model into the 
frequency dcmain. The next step is the formulation of the likelihood function 
of innovations and its maximization with respect to the unknown parameters. 

This step leads to the iterative scheme for parameter estimation, which updates 
the previous estimates by employing the second- and first-order gradients of 
the log-likelihood function. 

The linear airplane equations of motion are assumed in discrete-time form 

to be 


X(t+1) 

= $ x(t) + r u(t) + r w w(t) 

rt 

II 

o 

• 

• ., N - 1) 

(1) 

y(t) = 

H x ( t) + v (t) 

• 

o 

n 

-M 

. ., N - 1) 

(2) 
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E{x(0) } = 0 
E{x(0) x T (0) } = P 0 


( 3 ) 


where x(t) is a state vector, u(t) is a control vector, w(t) is a 
process-noise vector, y(t) is an output vector, and v(t) is a measurement- 
noise vector. 

It is assumed that 


(a) $, T, r w , and H are constant matrices 

(b) $ is stable 

(c) ($,T) and ($,r w ) are controllable pairs 

(d) ($,H) is observable 

(e) w and v are stationary, Gaussian uncorrelated noise sequences with 


E{w(t)} = 0 
E{w(t) w t (T)} = 


Q 6 t,T 


E{v(t)} = 0 
E{v( t) v T (X)} = 



E{v(t) w t (t)} = 0 for all t and T 


(4) 


(5) 


(f) N is even 

In the general case the unknown parameters will occur in the matrices 
r, r w , H , Q f R, x(0) , and Pq. Their estimation may be extremely diffi- 
cult because of the algorithm complexity (see ref. 15) and possible identifi- 
ability problems (see ref. 16). The system parameter estimation will be 
simplified by formulating a steady-state Kalman-f ilter representation of 
equations (1) and (2) and by considering the unknown parameters in this 
representation. 

The conditional expected value of the state vector is defined as 


x ( t) = E{x(t) |y(0) y (1 ) . . . y(N-l)} 


(6) 
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The innovations are defined as 


V(t) = y(t) - H x(t) (7) 

and the covariance matrix of state variables is defined as 

P = e{ Tx(t) - x ( t) ] l"x ( t) - x(t)] T > (8) 

Then the steady-state Kalman-f ilter representation of the system described by 
equations (1) and (2) is 

x(t+l) = # x ( t) + r u ( t) + K V(t) (9) 

y(t) = h x ( t) + v(t) (10) 


Reference 17 shows that the innovations V(t) form a sequence of independent 
Gaussian vectors with 


e{v ( t) } =0 1 


E{V(t) vT(T)} = B5 t , T J 

(ID 

The definition of the gain matrix K 
ence 15 or 16 in the form 

in equation (9) can be found in refer- 

K = $PH t B -1 

(12) 

B = HPH t + R 

(13) 

p = $p$t - kbk t + 

(14) 


For the further development of the identification algorithm all time func- 
tions in equations (9) and (10) are written in terms of their Fourier series 
expansions. As stated in appendix A, the Fourier series expansion of random 
variables holds in the mean-square sense. If the Fourier series component of 
x(t) is defined as 
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-I N-l 

x(n) x (t) exp(-jno) 0 t) (15) 

N t=0 


for 


N N N 

n ~ — If • • ■ , 0,1 , * • • f 1 

2 2 2 


2it 

where u)q = — , and similarly for the other variables, equations (9) and (10) 

N 

transformed from the time to the frequency domain have the form 

z n x(n) = $ x(n) + T u(n) + K \5(n) (16) 

y (n) = H x (n) + \5(n) (17) 

In equation (16) z n = exp(jnu)Q) , which follows from the relationship 
! N-l 

— ^ x(t+l) expt-jnojQt) = x(n) exp(jnu)Q) 

N t=0 


assuming that x(0) = x(N) = 0. As proved in appendix A, the transformed inno- 
vations \5(n) are uncorrelated, orthogonal, and Gaussian random variables with 


E{\5 (n) } = 0 


E{\5(n) V* (n) } 


>w 


N 


(18) 


where Sw is the spectral density of v(t), and v*(n) is the complex conju- 
gate of $ T (n). It follows from equations (16) and (17) that 

y (n) = H(z n I - $) -1 r u(n) + [H (z n I - $) _1 K + i] v (n) 

= T-j (n,©) u(n) + T 2 (n,0) V(n) (19) 
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where I is the identity matrix, T-) and t 2 are the system transfer func- 
tions defined as 


T-, (n,0) = H(z n I - $)-T (20) 

T 2 (n,0) = H (z n I - _1 K + I (21) 


and © is the vector of unknown parameters in equations (9) and (10). Equa- 
tion (19) is invertible in the sense that y(n) can be solved for directly 
in terms of $(n) and (n) in terms of y(n) (see ref. 18). This implies 
that T 2 is nonsingular. Therefore, from equation (19) 


\5(n) = T^ 1 y(n) - TjfTi u(n) 


( 22 ) 


To obtain the likelihood function, i.e., the joint probability density of 
the transformed innovations \5(n) (assuming that all parameters are known), 
a vector \5 m consisting of all innovations up to and including frequency m 
is introduced. Therefore 


\ - 





(23) 


Assuming that the probability distribution of \5 m has a density p[$ m ] , then 
it follows from the definition of conditional probabilities that 


pNJ -pRwlVil ptVi' 


(24) 


Repeated use of this formula gives the expression for the likelihood function 
as 


pr\] = pf^WlVll P^(m-l)|v m _2] 




(25) 


Because the distribution of \*(m) is Gaussian, then the distribution of $(m) 
given \>(m-l) is also Gaussian; i.e.. 
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exp[-N V* (m) Svv ^(m)] 

p[\>(m) |v(m)-l)] = 

TT r | Svv/N| 


( 26 ) 


as follows from the definition of a complex multivariate distribution in appen- 
dix B. In equation (26) r is the dimension of the innovation vector. 

Using equation (26) the logarithm of equation (25) can be written as 


N 

--1 
2 

J(0) = -N > \5* (n) S V v \5(n) - N log |Svvl “ Constant (27) 

N 

n=- - 
2 


In the log-likelihood function given by equation (27) , the unknown parameters 
are the elements of the matrices $ , T, H f K, and Sw An estimate of the 
unknown parameters is obtained by minimizing the log-likelihood function from 
the feasible set of parameter values. Optimizing the log-likelihood function 
for parameters in S^v gives 

Sw = z n V (n) \5*(n) (28) 


N 

—1 
2 

where £ n = The estimates of the remaining unknown parameters are given 

N 

n= — 

2 

by the root of the equation 


3j (0) 


30 


0 =© 


0 


(29) 


for Svv replaced by Sw This root can be found by a modified Newton-Raphson 
iteration (e.g., ref. 19) as 


0 = 0g + A0 


(30) 


n 



where the step size 


A© for parameter estimates is given by 


A0 



3J (0) 


3© 




( 31 ) 


The index 0 at the matrix M indicates that its elements were computed for 
0 = 0Q • In equation (31) M is the Fisher information matrix 


M = -E < 


3 2 J (0) 
30 30 T 


(32) 


Because the step size for parameter estimates is a vector with real elements 
only, and the log-likelihood function is real, the expressions for the first- 
and second-order gradients of J (0) are also real; i.e.. 


3J (0) 

30k 


-2N Re Z n V* (n) 



3v(n) 


(33) 


and 


3 2 J (0) 3v*(n) _i 3\>(n) 

= -2N Re E n S V v 

3© k 3©^ 3©£ 30 k 


(34) 


The expressions for the elements of the information matrix and the gradi- 
ent of the log-likelihood function are developed in appendix C in the form 


r_ * 
3 t 


MkS, = 2 Re Z n /Tr 


Ui l / *\-l _t UA 1 

i^r 2 ) SwT2 a5: Suu(n) 


3Ti 


30^ V 


+ Tr 


9 t 2 / *\-1 3 T 2 


( T2 ) 

3©n V Z / 




; 30£ 

1 

L 


To 

30 k 


(35) 


and 


3 J (0) 

3©k 


= -2 Re £ n <v Tr| 


9T 1 / * \- 
30b- V 2 ) 


„-l „-li 

9©k' 


Svv ^2 (Syy(n) ” T*j Syu (n) j — Tr| 


9t 


2 / * \-l 


if * v 
t 2 j 

3© k \ / 


( 36 ) 
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where 


S uu (n) 

A 

S yu (n) 


= N u(n) 
= N y (n) 


u* (n) 


u*(n) 


(37) 


are the estimates of input spectral densities and cross-spectral densities, 
respectively. 

The final estimates of unknown parameters have the following properties 
(see refs. 20 and 21): 

They are consistent; i.e.. 


lim P{ | 0 — 0| S e} = 1 

N-m» 


(with e arbitrarily small) 

They are asymptotically unbiased; i.e., 

lim E{0} = 0 

tf*oo 


and they are asymptotically efficient with 


E{ (0 -©)(©- 0) T } ^ -E 


3 2 J (©) 
30 3©T 


(38) 


Because of equations (32) and (38) , the inverse of the information matrix pro- 
vides the Cramer-Rao lower bounds on the variance and covariance of errors in 
the estimated parameters. 


OUTPUT ERROR METHOD 

If the process noise is zero and the initial states are assumed to be 
equal identically to zero, i.e., w(t) = 0, x(0) = 0, and P(0) = 0, the 

state-covariance matrix is also zero. Then, as follows from equations (12) 
and (21) , the Kalman gains are zero and T 2 = I. The innovations are reduced 
to output errors 


13 



v(n) = y(n) - T(n,©o) G(n) 


(39) 


where T(n,0Q) is equal to T-| (n,0) defined by equation (20). For 0(j 0 

the innovations v(n) -*• v(n) and S vv -»■ S w . The expressions for the 
elements of the information matrix and gradient of the log-likelihood function 
are obtained by simplifying equations (35) and (36) as 


and 


M kS , = 2 


Re 


N 

--1 

2 

SI 

N 

n=“- 

2 




3T 


^UU 


3j (©) 

~^T 


-2 Re Z n Tr| 


9t 

3®k 


s vvl s yu( n ) 


T s uu 



(40) 


(41) 


The expressions for the information matrix and the gradient of the log- 
likelihood function can also be easily derived from the simplified log- 
likelihood function, which takes the form of the output error cost function 


J (0) = -Nl n \5*(n) V(n) (42) 

These expressions are 

M(0) = 2N Re £ n A* (n) A(n) (43) 


3J (0) 
30 


-2N Re E n 


A* (n) S 


-1 

w 


V(n) 


(44) 


where A(n) is the sensitivity matrix whose elements are equal to 
3 [t ( n,©o) u(n)]/3©£. 

In some experiments airplane transfer functions are measured directly 
using a harmonic input or are determined from measured input-output time 
histories. Then the cost function includes a transfer function error rather 
than an output error. The cost function is therefore formed as 

J (©) = -N£ n [T E (n) - T(n,© 0 )]* S~J !T E (n) - T(n,© 0 )] 


14 


(45) 



where T is a vector which includes system transfer functions as elements. 
These transfer functions are computed from equation (20) for a given 0g. 

Both cost functions (42) and (45) can be minimized with respect to unknown 
parameters in $, G, and H or with respect to transfer function coefficients 
in T. The estimates are obtained from equations (31 ) , (43) , and (44) ; the 
spectral densities are given by equation (28) using pertinent residuals. 

For a system with a single input, the output error cost function with mea- 
sured transfer functions (frequency response curves) is defined as 

J(0) = -NZ n {u(n) [t e ( n) - T(n,© 0 )]}* S ^ 5(n) [t e ( n) - T(n,© 0 )] (46) 

In this formulation the scalar variable u(n) may be interpreted as a weight- 
ing function expressing the reliability of the measured data according to the 
harmonic content of an input. 


DISCUSSION 

The frequency domain identification has several features which are distinct 
from the time domain approach. They are mainly associated with the model repre- 
sentation and estimation algorithm. There is, however, the equivalence in the 
cost function used in the time and frequency domains as expressed by Parceval's 
theorem. This theorem postulates the relationship between the squared magni- 
tudes of the Fourier transform pairs. It therefore states that the time domain 
cost function. 


J TO = I t vT (t ) svi V(t) 


(47) 


N-l 

where Z t = , is equal to the frequency domain cost function, 

t=0 



Z n \>*(n) S v v ^(n) 


(48) 


Using equation (15) the frequency domain cost function can be written as 


1 _i 

J FD = “ Z n S t vT (t) exp(jnu) 0 t) S V v V ( T ) exp(-jnO) 0 T) 


1 


= - £ t Z x ^ T (t) Svv V ( T ) £ n exp T jnug ( t - T) ] 
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N-l 

where Z T = . But according to appendix A, 

x=0 


£ n exp[ jnu>Q (t - X)] = N 

(for t - X) 

= 0 

(for t/x) 


Therefore, 


JpD = ^t vT (^) 


8 $ 


V (t) = jrpj) 


The equivalence of both approaches is no longer valid if the frequency domain 
cost function is restricted to a given frequency range. Such a restriction is 
not necessary, but it is an option which is a strong point in favor of fre- 
quency domain analysis with respect to time domain analysis. The selected 
frequency range of interest was used, for example, in reference 9, where air- 
plane rigid modes were separated from elastic ones. For similar results in the 
time domain the data must be filtered accordingly. 

The early airplane estimation techniques in the frequency domain were 
using measured frequency response curves only. This approach could have an 
advantage when repeated measurements under the same conditions are available. 

A hypothesis concerning the model adequacy can be tested using the variance 
estimates from scatter around the mean and from residuals (ref. 5) . On the 
other hand the simultaneous analysis of repeated maneuvers for obtaining a 
single set of estimates with increased accuracy can also be applied to directly 
measured or transformed time histories, in general, transformed input-output 
time histories are preferred in frequency domain parameter estimation. The 
inaccuracies of frequency response curves computed from transformed inputs and 
outputs can be quite pronounced for frequencies in which the harmonic content 
of an input is close to zero. 

The transformation of model equations into the frequency domain replaces 
differentiation and convolution with multiplication. As a result the sensitiv- 
ity equations in the nonlinear estimation algorithm are reduced to uncoupled 
algebraic expressions. This simplification can be appreciated mainly in cases 
for which convolution integrals are included in the equations of motion 
(ref. 11). 

The computational differences between the time and frequency domains dis- 
cussed so far could be viewed as advantages of the frequency domain analysis. 
There is, however, a substantial disadvantage of the airplane identification 
in the frequency domain. This approach is limited, for practical reasons, to 
only linear equations of motion with constant coefficients. The computing time 
needed for parameter estimation in the frequency domain (transformation of mea- 
sured data included) is about 50 percent more than in the time domain. The 
assessment was obtained from the number of equations used in both domains for 
one iteration when the algorithms were applied to the system of equations with- 
out convolutions and process noise. 
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The estimation algorithm was developed for a linear discrete-time model. 
The airplane equations of motion are, however, usually given in a continuous 
form as 


x = Fx + Gu + G w w 


(49) 


where the unknown parameters can be in the matrices P, G, and G w . For the 
continuous model (eq. (49)), the expressions for the information matrix and 
gradient of the log-likelihood function remain the same as equations (35) 
and (36) . But now the transfer functions are defined as 


T-, ( 03 , 9 ) = H(jU)I - F) _1 G (50) 

T 2 ((O,0) = H( jtOI - F) -1 K C + I (51) 


where the Kalman-f ilter-gain matrix is obtained from the relationships (see, 
e.g., ref. 22) 


K c = PH t R _1 


and 


FP + pfT _ ph t R - 1 HP + (^QcG^ = 0 


In the model formulation it was assumed that the initial conditions were 
equal to zero, that the model described a stable motion of an airplane, that 
there were no a priori known values of stability and control parameters, and 
that the measurement noise was Gaussian and uncorrelated. If the initial con- 
ditions differ from zero, the additional term H(z n I - x(0) would have to 

be included in equation (19) or the new transfer function H(jo)i - F)~^ x(0) 
would have to be added to those defined by equations (50) and (51). Then 
the vector of unknown parameters can be augmented by the vector of initial 
conditions. 

If the airplane motion includes an unstable mode, the parameter estimation 
still can proceed provided that the degree of instability is not high. A large 
instability, on the other hand, can result in excessive transient motion due to 
nonzero initial conditions and/or the input and thus limit the validity of the 
linear equations of motion. If the a priori mean values and variances of some 
parameters are known, they can be included in the estimation procedure. In 
this case the cost function must be expanded in a similar way as indicated in 
reference 19. 

The maximum likelihood method developed earlier assumed a Gaussian, uncor- 
related measurement noise. If the random sequence representing this noise is 
correlated, the estimation algorithm does not change. The constant values of 
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spectral densities S vv or S vv are merely replaced by the frequency depen- 
dent values estimated from expressions similar to equation (37) . 


EXAMPLES 

As examples the parameters of a small general aviation airplane were esti- 
mated from computer-generated data and from measured data in still and turbu- 
lent air. (Some of the data for examples 1,3, and 4 are from ref. 14.) For 
all examples the model of the airplane was based on simplified longitudinal 
equations of motion with the atmospheric turbulence (gusts) approximated as a 
Gauss-Mar koff process of first order. The model equations (continuous form) 
are developed in appendix D. When the state and output equations (D4) and (D9) 
are transformed into the frequency domain and rearranged, they have the form 


xD = G<5 e + G w w 

(52) 

y = Hx + v 

(53) 


The state and output vectors are specified as 


x = [a q Wg] T 
y = [a v q a z l T 


where the random input is 
with E{w} = 0 and E{w 2 } 
formulated as 

assumed to be a Gaussian, uncorrelated noise process 
= The matrices in equations (52) and (53) are 



'jo) - k-,C Za 

-(l + k 2%) 

" k 3 C Z a 



D = 


3 “ ~ k 6 c mq 





0 

0 

jd) + c-| 



G = 

_ kl %e "^e °] T 



r -T 

G w = 0 -k 8 c 2 C§ g c 2 
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kg 


- k io 


kg 

V 0 


v 0 v 0 

— jo) 0 

g g 


The known constants k-j , k 2 / • . . , k-jg^ c-j , and C 2 , and the pitching- 

moment derivatives C^, C^, C ^q' and ^5 are def ^ ned in appendix D. 

Example 1 

The output data were computed from equations (D4) and (D9) for given 
inputs 6 e and w (with o g = 1 m/sec) and for a given set of parameters. 
To the computed time histories of output variables, an uncorrelated and 
Gaussian measurement noise was added. The measurement-noise standard errors 
were selected as 


a a = 0.0028 rad o Q = 0.0063 rad/sec a a = 0.02g 

^ z 


The time histories of the input and output variables are plotted in figure 1 . 
For the estimation algorithm these data were transformed into the frequency 
domain using the Filon integration formula. The sampling interval of the 
transformed data was 1.047 rad/sec. The transformed data were truncated at 
the frequency interval ±20.944 rad/sec because outside this interval their 
amplitudes were very small. 

The steady-state Kalman filter representation of the airplane motion 
described by equations (52) and (53) is 


xD = GS e + K C V (54) 

y = Hx + V (55) 


It was assumed that the parameters c-j and C 2 * the initial conditions, and 
the variances of the process and measurement noise were known. Also assumed 
as known was the parameter because of the identif iability problem. 
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The last assumption is substantiated by the small effect of the term kgo-jCg 
on the airplane motion. The vector of unknown parameters was therefore ^ 
formed as 

0 = [ C Zct C Z 5e C Zq ^ ^ C m 6e K H K 12 • • • ^ 33 ] 


where Kn , K 12' • . .» and K33 are the elements of the Kalman-fi lter-gain 

matrix K c . The initial values of these elements were computed from equa- 
tions (50) and (51) . 

First, the unknown parameters were determined by the maximum likelihood 
method developed. In table I the estimated stability and control parameters 
are compared with their true values, and the estimated Kalman gains are com- 
pared with their initial values. The agreement between the first set of param- 
eters is, in general, very good. The estimated values of the Kalman gains 
differ significantly from their initial values, and the standard errors (lower 
bounds) of these parameters are quite high. This indicates low accuracy of 
these estimates. When, however, the Kalman gains were fixed on their initial 
value, the estimates of the airplane parameters were farther from the true 
value, as indicated by results in the fourth column of table I. The last set 
of airplane parameters was obtained by considering no process noise effect on 
the output data. These estimates are also less accurate than those obtained 
by the maximum likelihood method with all 15 unknown parameters. Table I also 
includes the variance estimates of the residuals. The limited experience 
obtained from this example indicates that for stability and control parameter 
estimation from data with pronounced effect of the process noise (i.e., 

Og £ 1 m/sec), the algorithm in its complete form should be used and the Kalman 
filter gains should be treated as the additional unknown parameters. 


Example 2 

In this example the measured data in turbulent air were used in the same 
model as in the previous example. The measured input-output time histories are 
presented in figure 2. For the parameter estimation the transformed data were 
taken from the frequency interval ±9.817 rad/sec. The standard error of the 
vertical gust velocity was determined from the part of the measured data with 
<5 e = Constant to be c?g = 1.12 m/sec. The values for measurement-noise stan- 
dard errors were taken Iran the results in reference 23 as 


a a v = 0.0017 rad = 0.005 rad/sec = O.Olg 

The estimated stability and control parameters are given in the third and 
fourth columns of table II. In the first case (fourth column) the Kalman gains 
were treated as unknown parameters; in the second case (third column) they were 
set equal to zero (assumption of no process noise) . The inclusion of the pro- 
cess noise in the model resulted in better accuracy of the parameter C^, as 
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indicated by its comparison with the average value obtained from the estimates 
in the time domain (see ref. 23). On the other hand, the process-noise consid- 
eration in the estimation process degraded the estimates of the parameters C z 
and Cm. No explanation for this degradation could be found. q 


Example 3 

From the measured time histories in still air which are presented in fig- 
ure 3, the transformed input and output data and the frequency response curves 
relating all three outputs to the elevator deflection were obtained. By setting 
w » 0 the state vector in equations (52) and (53) was changed to x = [ot,q] T 
and the matrices D, H, and G were simplified accordingly. The unknown 
parameters were estimated from the minimization of the cost function, given by 
equation (42) for the transformed data and by equation (45) for the frequency 
response curves. 

The estimated parameters are given in the sixth and seventh columns of 
table II, and they are compared with the results from the time domain estima- 
tion given in the fifth column of the same table. The three sets of estimates 
from the same flight agree well. The standard errors of the estimates in the 
frequency domain are, however, higher than those in the time domain. This 
could be due to truncation of the transformed data and additional inaccuracies 
in measured frequency response curves caused by taking the ratios of two com- 
plex numbers. The transformed data and those computed are plotted in figure 4; 
the measured and computed frequency response curves are plotted in figure 5. 

Both figures indicate some modeling errors in the equation for a v . It is also 
apparent from figure 5 that the measured frequency response curves are inaccu- 
rate around the frequency 6.4 rad/sec as a result of the low harmonic content 
of the input at the same frequency. 


Example 4 

The response of the airplane to turbulence was measured in two flights 
(designated run 1 and run 2 in table II) with the minimum pilot interference 
(6 e * 0) • From the time histories of the measured output variables the spec- 
tral density of the vertical gust velocity S w w and the cross-spectral den- 

sities S w ^q and s w g a z were computed. They are related by the airplane 

transfer function resulting from equations (52) and (53) . 

The state and output equations were modified in the following way: 

(a) In equation (52) w and S e were set equal to zero 

(b) Wg was assumed as a known input 

(c) The term kgc-iCg^Wg was replaced by -jookgCg^Wg 

(d) In equation (53) was set equal to zero 
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The matrices D, G, and H were therefore changed as 


jo> - k-,C Za -(l + k 2 c Z q ^ 


G = 




v 0 v 0 

— j(0 

g g 


and the model was formulated as 


S w GSy* M 

W 9' X ‘ w g w g 

s w g y = + ▼ 

where 

s w g x = ^ s w g a s w g qJ 


S w g y 



w g a z 


] 


T 


( 56 ) 

( 57 ) 


The vector of unknown parameters in these equations is formed as 


© = 




The estimated values of the first four parameters are given in the last two 
columns of table II. From the estimates of Cj|^ and the value of 
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was computed from equations (D6) and (D7) and included among the unknown param- 
eters. The agreement between the results from both runs is very good. The 
parameters also agree with the estimates from the still air measurement with 
the exception of the parameter c z a * This parameter has a smaller value than 
expected, probably because of some modeling errors in equations (56) and (57) . 

The measured spectral and cross-spectral densities from run 1 and those 
computed by using the estimated parameters are plotted in figure 6. It was 
verified that the large fit error in the phase of the cross-spectrum S w q 
did not affect the values of the estimated parameters significantly. Th§ esti- 
mate from turbulence and measurement demonstrates a possibility for using these 
data also for airplane stability and control parameter estimation. For certain 
model formulations the derivative of pitching-moment coefficient with respect 
to the rate of change in angle of attack can be estimated explicitly. 


CONCLUDING REMARKS 

A frequency domain maximum likelihood method has been developed for the 
estimation of airplane parameters from measured flight data. A discrete-type 
steady-state Kalman filter was used in the derivation of the computing algo- 
rithm. The time variables in the model equations were transformed into the 
frequency domain by using a Fourier series expansion. If the initial data 
were Gaussian and uncorrelated, the transformed data formed a complex random 
sequence which was uncorrelated, orthogonal, and Gaussian. Then, the likeli- 
hood function could be formulated as a multivariate distribution of complex 
innovations. 

The connection between the continuous form of airplane equations of motion 
and the developed algorithm is easily established. The algorithm can be sim- 
plified to the output error method with the measured data in the form of trans- 
formed time histories, frequency response curves, or spectral and cross-spectral 
densities. In general, transformed input-output time histories are preferred 
in frequency dcmain estimation. The inaccuracies of frequency response curves 
computed from transformed inputs and outputs can be quite pronounced for fre- 
quencies in which the harmonic content of an input is close to zero. The fre- 
quency domain approach simplifies the estimation procedure by reducing the sen- 
sitivity equations to simple algebraic expressions. It also provides an easier 
way than the time domain for using the data within a frequency range of inter- 
est. The serious disadvantage of the frequency domain identification is in its 
practical limitation to a system described by linear equations of motion with 
constant coefficients. It was shown that the cost functions in time domain and 
frequency domain approaches are equivalent. It is therefore necessary to con- 
sider both approaches as complementary and not contradictory. 

The maximum likelihood method has been applied to computer-generated and 
real flight data for the longitudinal motion of a small general aviation air- 
plane. In the first case the estimates obtained were more accurate than those 
from the simplified output error method, which did not consider the effect of 
the process noise. In the second case the results were inconclusive because of 
an insufficient amount of measured data. Then, the simplified algorithm was 
used with the flight data from measurements in still air and turbulent air with 
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no pilot input. The first set of results for deterministic input showed the 
expected similarity in parameter values obtained from the time and frequency 
domains. The estimates from turbulence measurements demonstrated a possibility 
for using these data also for airplane parameter estimation and for explicit 
estimation of the pitching-moment derivative with respect to the rate of change 
in angle of attack. 


Langley Research Center 

National Aeronautics and Space Administration 
Hampton, VA 23665 
April 1,1 980 
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APPENDIX A 


FOURIER TRANSFORM OF A STOCHASTIC PROCESS 
Let z ( t) , where t = 0 , 1, . . . , and N - 1 , be a real random sequence 

with 


E{z ( t) } = 0 


E{ z (t) z (T) } = R zz (t-T) 


(Al) 


Theorem 1 : If z(t) is periodic in the mean-square sense, then it can be 

expanded into a Fourier series: 


N 

— 1 
2 

z(t) = z(n) 

N 

n=-- 

2 


e jntogt 


(A2) 


where u)g = 2tt/N; and the coefficients z(n), given as 


z(n) 


i N-l . 

- z(t) e - ‘ ]nU ^ t 

N t=0 


(A3) 


are uncorrelated and orthogonal randan variables such that 

~N 

E{ z (n) } = 0 

S zz (n) > 

E{z(m) z* (n) } = — 6 m , n 


(A4) 


where S zz (n) is the spectral density of z(t) . 

Proof ; From equations (Al) and (A3) the expected value of z(n) is found as 


E{z (n) } = - E{z (t) } e" jnt ° ot = 0 

N t=0 
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To prove the second part of equation (A4) , the conjugate of z(n) is first 

multiplied by z(t) and then the expected value is taken; i.e.. 


E{z*(n) z(t)} 



j N-l 

E< > z(s) 
) s=0 



-EE{.(t) z(s) } e jna) 0 s 
N s=0 


1 N-l 

= -ZlRzz(t-s) 

N s=0 


e jna) 0 s 


The new variable T = t - s and the relationships 


N 

--1 
2 


Rzz CO 


= 1 Szz(n) e^ nW 0 T 


(A5) 


N N 
n=-- 


N-l 


S zz (n) -^ZZRzzCO e _3nW ° T 


(A6) 


T— (N-l ) 


are introduced (see, e.g., ref. 6). Then 


E{z(n) z ( t) } = - R zz (0 e jnW ° (t T) 

N t=-(N-l) 


1 e j na) 0 fc 

N 


S 2 z ( n ) 


(A7) 


Using equations (A3) and (A7) , the expected value of z(m) z*(n) can be 
formulated as 


•, N-l 

(m) z* (n) } = - EZ 
N t=0 


E{z 


S, z N-l 

* (n) z(t) } e “3-“0t 

m 2 t=0 


N-l 

e ~D (m-n) w 0 t 

N 2 t=C 


(A8) 


26 


APPENDIX A 


For m = n there is 

y 1 e -j(m-n)o) 0 t = N 
t=0 

For m ^ n the difference m - n = a, where a is an integer; therefore. 


-jaU)Qt 1 - e **** ^ - cos 2aTT + j sin 2air 

^0 6 ! . e ’ jWot 1 - e ' ja) 0 t 

and the proof of equation (A4) is thus completed. 

To prove equation (A2) , it is sufficient to show that the sequence z(t) 
tends to £ n z(n) exp(jnU)gt) in the mean-square sense; i.e.. 




E< 


z(t) - 


N 

— 1 
2 

z: 

N 

n= — 
2 


z (n) e 


jno) 0 t 


>-• 




(A9) 


The square in equation (A9) can be written as a product of itself by its 
conjugate. Equation (A9) is therefore changed to 


E{|z(t)| 2 } - Z n E{ z* (n) z ( t ) } e ^ nU) 0 t - £ n E{z(n) z*(t)} e ^ nC ° ot 
+ £ n Z m e{ z (m) z*(n)} e^ m n ^ a) 0 t _ o 


In the double summation all the terms with m ^ n are equal to zero. Using 
equations (A! ) , (A2) , and (A7) , the previous equation can be expressed as 


^zz (0) 


l S n S zz (n) 
N 


e jnw 0 t e -jna) 0 t 


R ZZ (0> + l S zz (n) 
N 


o 


Each sum above equals R zz (0) (see eq. (A5)); hence the whole expression 
equals Z ero and equation (A9) is proved. 
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Theorem 2 : Let z(t), t = 0, . . . , and N - 1 , be mutually stochastically 

independent random variables having, respectively, Gaussian distributions with 
e{z ( t) } = 0 and E{z 2 (t)} = a 2 . Then the sequence 


■. N-1 . 

z (n) = 1 z(t) e -3nw 0t 

N t=0 


(A3) 


where 


n 


N N 

■■ f • • • / 0 f 1 f • • • / “ 1 

2 2 


consisting of real parts z R (n) and imaginary parts z 1 (n) is a complex 
Gaussian and uncorrelated sequence with 


and 


E{z R (n) } = Eiz 1 (n) } = 0 


(AT 0) 


r 

i 

£ 

IN 

; i 

z R (n) 

T^ 


s zz 

0 

\ 

z 1 (m) 

z I (n) 

-J 

2N 

0 

s zz 


J m,n 


(All) 


Proof : It has already been proved that 

E{z (n) } = 0 


This result combined with the definition of the expected value of a complex 
random variable (eq. (B1)) implies that 

E{z R (n) } = E{z I (n) } = 0 

It has also been shown that (in eq. (A8)) 


E{z(m) z (n) 


S 77 (n) N-1 . . _ 

} = z e -j (m-n)a) 0 t _ 


>zz 


(n) 


N^ 


t=0 


N 


J m,n 


(AT 2) 
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Using the same approach as that for the development of equation (A8) 
be shown that 

E{z(m) z*(n)} = E{z R (m) z R (n)} + E{z J (m) z*(n)} 

+ jE{z I (m) z R (n)} - jE{z R (m) z I (n)} 

and 

E{z(m) z(n)} = E{z R (m) z R (n)} - Eiz-^m) z-^n)} 

+ jE{z R (m) z I (n)} + jE{z I (m) z R (n)l 
From equations (A8) and (A12) to (A14), two sets of equations can be 

E{z R (m) z R (n)} + Etz^fm) 

E{z R (m) z R (n)} - Eiz^fm) 

and 

E{z I (m) z R (n)} - E{z R (m) 

E{z R (m) z^n)} + Eiz^m) 

These two sets give the solution 

E{z R (m) z R (n)} = Eiz 1 (m) 

E{z R (m) z 1 (n) } = Etz 1 (m) 


z 1 (n) } = 


&zz ( n ) 




N 


J m,n 


z 1 (n) } = 0 





^zz 

z* (n) } = 6 

2N 


m,n 


z R (n) } = 0 


it can 

(AT 3) 

(A14) 
formed as 

(A15) 

(A16) 

(A17) 

(A18) 


which proves the validity of equation (All). Equations (A15) and (A16) 
have been developed for m and n being positive. In the case where 
(-N/2) ^ m f n ^ (N/2 - 1) the proof based on these equations is still valid 
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The only change might occur in equations (AT 5) , where the right-hand sides 
would be interchanged. The form of the solution given by equations (A17) 
and (A18) does not change. 

Finally, it is necessary to prove that z(n) is Gaussian. For this 
proof the concept of the moment-generating function is used. When z(n) 
is expressed as 

N-l 

z(n) = c t z(t) 
t=0 

where 


= I e -jn»„t 
C N 


then the moment-generating function of the variable c b z(t) is given accord- 
ing to reference 24 as 


E 


bctz (t) 


exd 


a 2 

-(bc t ) 


*(bc t ) 


(A19) 


where b is a constant independent of c t z(t). Thus, the moment-generating 
function of z(n) is 


E{exp bfco z(0) + c-] z(l) + . . . + c^-j z(N-l)]} 


= E{exp bCQ z(0)} E{exp bc-j z(l)} 


E{exp bcN-i z(N-l)} 


N-l 

a 2 * 


r b 2 a 2 * i 


b 2 a 2 

II exp 

— (bc t ) (bc t ) 

= exp 

c t c t 

= exp 


t=0 

2 


2 t=0 


2N 


which means that z(n) is Gaussian with zero mean and variance a 2 / N or, 
using equation (A4) , with variance S zz /N. The sequence z(n) is formed 
by a collection of uncorrelated Gaussian random variables. 
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COMPLEX RANDOM VARIABLE 

A complex random variable z is a complex number z(C) determined by 
an outcome ? ; i . e . , 


z(C) = z R (C) + jz I (C) 


such that the functions z R and z 1 are random variables. By definition, the 
expected value of a complex random variable is 


E{z} = E{z R } + jE{z 1 } (Bl) 

the variance is 

°z = b{|z-e{z}| 2 } (B2) 


and the covariance is 


e{ [z^-eCz^}] [z£-e{z$}]} 


If 


(B3) 


Elz^zJ} = eCz^} e{zJ} (for k / S,) 

then the complex random variables z-j , Z2, . . . , and z n are uncorrelated. 

They are orthogonal if 

E{z k z$} = 0 

If the complex random variable z = z R + jz 1 has its real and imaginary 
parts normally distributed with 

e{z r } = e{z 1 } = 0 


e{ (z r ) 2 } = e{ ( z 1 ) 2} = “ CT z 
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and these parts are stochastically independent, then the distribution of a 
complex random variable z will be defined as a joint distribution of the 
independent variables z R and z 1 such that 


p[z R ] plz 1 ] = exp 

” (z R ) 2 /o z 

1 

exp 

1 

<N N 

CM 

H 

N 

s * 

fi° z 


f°z 



expl-z 


tto. 



(B4) 


If z is a complex vector of dimension r, then the Gaussian distribution 
is defined as 


pfz] = exp(-z*Z -1 z) 

7T* \E 


(B5) 


where 2 is the covariance matrix of z. This is a Hermitian nonnegative 
square matrix of dimension r. If z has components that are stochastically 
independent, then 2 is a diagonal matrix. Equation (B5) agrees with the 
definition of the Gaussian distribution presented in references 24 and 25. 
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INFORMATION MATRIX AND GRADIENT OF LOG- LIKELIHOOD FUNCTION 

From equations (32) and (34) the elements of the information matrix are 
given as 


9v* _i 9v 


K kl = 2N Re E^Z n — S vv — 

90k "S, 


where the expected value can be written as 


E 








(Cl) 


(C2) 


The measured outputs and innovations are given by equations (19) and (22) as 


y(n) = T*j u(n) + T 2 V(n) (C3) 

and 

\5(n) = T 3 y (n) - T 4 u(n) (C4) 


where T 3 = T 2 ^ r T 4 = T^T-j, and T-j 
and (21 ) . Therefore, 


and T 2 are defined by equations (20) 


E 





= T 3k E{yy*} T^ - T 4|t E{uy*} T^ 

- T 3k E{yu*} T 4a + T 4k E{uu*> T 4& (C5) 
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where 


T 3 k 


9t 3 

3®k 



3t 


* 

4 




Because u(t) and v(t) are uncorrelated, the transformed variables u(n) 
and v (n) are also uncorrelated. Using equations (C3) and (B4) the expected 
values in (C5) can be expressed as 


'N 


E{yy*} = E jj-,uu* T* + = ~ (TiS uu T* + T 2 S VV T 2 ) 


E{uy*} = E{ uu*T* + u\>*T 2 ) = - 

N 

E{yu*} = E{T-,uu* + T 2 3 u*} = - T-j S uu 

N 


(C6) 


E{uu*} 



J 


After substituting equation (C6) into equation (C5) and some tedious 
manipulation, 


9v* 3v l i -1 * * _i 

E< 90^ 9%^ = T2 + * 2 T 1k S uu T U< T 2> 


(C7) 


From equations (C2) and (C7) , 


9v* 3v { 1 

E < S vv 

90 k W 905, ( N| 


Tr^T* JL (T 2 ) -1 S v i(T 2 ) _ ' 1 T-, k S uu j + Tr ^ (T 2 ) 


(C8) 
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Substituting equation (C8) into equation (Cl), the final form for the elements 
of the information matrix is obtained as 


— 2 Re Z 




~^ V W 



+ Tr 






(C9) 


From equation (33) the element of the gradient of the log-likelihood 
function is 


3J (©) 

a© k 


-2N Re Z n v*S 


_1 3v 

vv 


80 k 


= -2N Re Z 


n 


3v* 

30 k 


3 vv^ 


(CIO) 


Using equation (C4) 


3\5_ 

30 k 


t 2 l ?2 a 



T l k 3 


(Cll) 


Then 


3 \ 5 * 

3© k SvvV 


" ( T 2 k T 2 v) S v iv - (^T^uj S v Jv 


(Cl 2) 


where 


(T 2k T 2 v) S"iv = Tr|^VV*SviT$ k (T^)-lJ 


(Cl 3) 


and 


( T 2 lT l k a ) S viv - Tr£yu*T* k (T 2 ) - 1 S^T 2 ] J - Trj^uu*T-, k (T 2 ) _ ' I S^T 2 1 T-, J 


(Cl 4) 


35 


APPENDIX C 


Introducing 


* 1 

VV = — Sxjxj 

N 


* 1 - 
yti = - s yu 

N y 


* 1 ~ 
uu = - S uu 
N 


and approximating = I, equation (Cl 2) is changed as 


3\5 


* 


a® k 


^ = 1 




(T2)" 1 



Finally, substituting equation (Cl 5) into equation (CIO), the elements in 
gradient vector are obtained as 


3J(0) 

3©k 


-2 Re Z n < Tr| 


3t* , , A 

—— (T2)“ 1 S y v T 2 (Syu 
dW k 


TtSuu) 


Tr 


3 To 

^- (t 2 } 

3© k 


-1 



(Cl 5) 
the 

(C16) 
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EQUATIONS OF LONGITUDINAL MOTION OF AN AIRPLANE IN PRESENCE OF TURBULENCE 

The airplane equations of motion are referred to the body axes. They are 
perturbed equations for datum conditions corresponding to steady horizontal 
flight. The equations are based on the following assumptions: 

(1) The airplane is a rigid body 

(2) The elevator deflection and turbulence excite the longitudinal motion 

during which the airspeed remains constant 

(3) The turbulence is approximated by a one-dimensional gust field, and 

the angle-of-attack and pitch-rate perturbations due to turbulence 
are given as 



and 

q g = -Og 


(4) The aerodynamic model equations for the increments in C z and C m 
have the form 


Ac Z = C Za (a + Og) + C Zq (q - a g )^- + C Z6e 6 e 


c c 

Ac m - <**<“ + o g ) + + a g )— + C^q - a g )— + C ra5e 6 e 


where the input and output variables are the increments with respect to the 
initial steady flight. 
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Using these assumptions, the longitudinal equations of motion can be 
expressed as 


a g c 


i - 9 + ^1%“ + % S" + C ^e 5 e + c %“g - % ' 9 Sin 9 ° 9 <D1) 


a g c 


qSc / ac qc 

4 = + + + + Vg - % — + ^ 


a g c 


2V 0 2V 0 


(D2) 


where g sin 0 q 9 in equation (D1 ) is negligible. 

In equations (D1) and (D2) ctg is a stochastic variable. Its spectral 
density can be modeled, for example, by a Dryden formula (see ref. 26) . In the 
further development it will be assumed that the turbulence velocity component 
Wg is a Gauss-Markoff process of first order governed by the differential 
equation 


Wg = -c-)Wg + C2W (D3) 

where 

Vo 

Ci = 2.4 — 

L g 

c 2 = v 0 



L„ is the scale of turbulence, and w is the uncorrelated noise process with 
Elw) = 0 and e(w 2 } = Cf2. 

When equation (D1 ) is substituted into equation (D2) and equation (D3) is 
considered, then the state equations for the longitudinal motion of the air- 
plane will have the form 
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• 

a 


klCz bt 

( 1 + k 2 C Z q ) 

k 3 c ^ + k 4 c l c Zq 

a 


k l c Z5e 


-k 4 c 2 c Z q 

• 

q 

= 

ksC % 

k 6 c mq 

k 7^'m x + k 8 c l^mq 

q 

+ 

k 5 c n^ e 

5 e + 


_*g_ 


0 

0 

-c-, 

w g_ 


0 


c 2 


(D4) 


where 


psv 0 

k 1 = 

2m 

pSc 

k 2 - . 

4m 

„ ps 

k 3 = — 

2m 

pSc 

k 4 = 

4mV 0 

2- 

pSV 0 c 

pSV 0 c 


k 5 = 


2 1 y 


k 6 = 


4Iy 


k 7 


pSV 0 c 
2 Iy 


pSc 2 



(Note: 


k 4 c 2 c Z q * °) • 


. pSc 

" °% + 4m °% Cz a 


(D5) 


^ ~ C m q 
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(D7) 


, PSc 

= Cm 6e + ^"a Cz 6e 


(D8) 


For the parameter estimation from measured data, state equations (D4) are 
completed by the output equations 




(D9) 
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TABLE I.- ESTIMATED AIRPLANE PARAMETERS FROM COMPUTER-GENERATED DATA 


Parameter 


Initial value 


Estimate of parameter 
(a) 


With process noise 






no process noise 



K estimated 

K fixed 


C Z« 

-5.0 

-5.07 

-5.72 

-5.2 

LX 


(.091) 

(.20) 

(.10) 

C Z CT 

-20 

-23.8 

-17 

-15 

q 


(2.3) 

(2.7) 

(4.0) 


-1 .0 

-1.17 

-1 .1 

-.57 

06 


(.071) 

(.17) 

(.31) 


-.80 

-.82 

-.94 

-.83 


(.015) 

(.031) 

(.014) 

C^a 

-24 

1 

to 

>(>> 

• 

o 

-25.5 

-22.1 

iuq 


(.52) 

(.38) 

(.65) 


-3.3 

-3.21 

-3.16 

-3.18 

0 6 


(.062) 

(.080) 

(.050) 

K n 

-.0563 

-.44 

-.0563 

0 



(.oil) 



*12 

2.544 

2.31 

2.544 

0 



(.010) 



*13 

.0025 

-.21 

.00 25 

0 



(.010) 



*21 

2.438 

1 .84 

2.438 

0 



(.044) 



k 22 

9.7 

10.8 

9.7 

0 



(.14) 



*23 

-.018 

.09 

-.018 

0 



(.97) 



k 31 

213.0 

156 

213.0 

0 



(8.9) 



k 32 

-5.065 

-8 

-5.065 

0 



(2.1) 



*33 

-12.63 

-15 

-12.63 

0 

2 


(6.2) 



S$ a X 10 5 

— 

.67 

.72 

.70 

2 





aci x io 4 

V q 

— 

.30 

.60 

.82 

2 





9 \) a 1 
a z 


.50 

1 .1 

1 .6 


Assumed 


a Numbers in parentheses are Cramer-Rao lower bounds on standard errors. 
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TABLE II.- ESTIMATED AIRPLANE PARAMETERS FROM MEASUREMENTS IN TURBULENT AND STILL AIR 


Estimate of parameter 
(a) 



umbers in parentheses are Cramer-Rao lower bounds on standard errors, 
ran the maximum likelihood estimates in the time domain (ref. 23 ) . 
ransformed input and output time histories. 

^Frequency response curves. 
e Spectral and cross-spectral densities. 

^Computed fran estimated and C^. 
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Figure 4.- Concluded. 
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Figure 5.- Continued. 




20 


15 


a 

z 



,10 


I e I 

g units/ rad 


5 


0 


-45 


<P 


a 6 ' 
z e 

deg 


-90 


-135 

-180 


54 





Figure 6.- Continued. 
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Figure 6.- Concluded. 
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